{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook Set C\n",
    "\n",
    "---\n",
    "\n",
    "## 04 - Investigating the DFT\n",
    "In this notebook, we investigate the DFT by first undertaking a simple example and then exploring an variety of DFT topics. These include DFT periodicity and symmetry, positive and negative frequencies, frequency representation, spectral leakage, windowing, scalloping loss, and zero padding.\n",
    "\n",
    "## Table of Contents\n",
    "* [1. Introduction](#introduction)\n",
    "* [2. An Example of the DFT](#an_example_of_the_dft)\n",
    "* [3. Magnitude and Phase](#magnitude_and_phase)\n",
    "* [4. Plotting the DFT](#plotting_the_dft)\n",
    "* [5. Positive and Negative Frequencies](#symmetry_and_periodicity)\n",
    "* [6. Magnitude, Power, and the Logarithmic Scale](#magnitude_power_log)\n",
    "* [7. Frequency Bins and Spectral Resolution](#frequency_bins)\n",
    "* [8. Spectral Leakage](#spectral_leakage)\n",
    "* [9. Windowing](#windowing)\n",
    "* [10. Scalloping Loss](#scalloping_loss)\n",
    "* [11. Zero Padding](#zero_padding)\n",
    "* [12. Conclusion](#conclusion)\n",
    "\n",
    "## References\n",
    "* [1] - [f. j. harris, \"On the use of windows for harmonic analysis with the discrete Fourier transform,\" Proceedings of the IEEE, Vol. 66, No. 1, Jan. 1978, pp. 51-83.](https://doi.org/10.1109/PROC.1978.10837)\n",
    "\n",
    "## Revision\n",
    "* **v1.0** | 20/01/23 | *First Revision*\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "We will begin by creating a discrete waveform that we can use as test stimulus for the DFT. Below, we define the regular sampling period $t_s$ and the number of discrete samples $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "fs = 48e3 # Sampling frequency\n",
    "ts = 1/fs # Regular sampling period\n",
    "N = 16    # Number of samples\n",
    "n = np.arange(N) # Sample points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create two cosine waves that we will add together to produce the final test waveform.\n",
    "\n",
    "Waveform $a$ has a frequency of 3kHz and waveform $b$ has a frequency of 9kHz.\n",
    "\n",
    "$$\n",
    "a(n) = \\frac{1}{2}\\cos(2\\pi 3000 n t_s + \\pi/4) \n",
    "$$\n",
    "\n",
    "$$\n",
    "b(n) = \\cos(2\\pi 9000 n t_s)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0.5*np.cos(2*np.pi*3e3*n*ts + np.pi/4)\n",
    "b = np.cos(2*np.pi*9e3*n*ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding each cosine wave together produces $x(n)$.\n",
    "\n",
    "$$\n",
    "x(n) = \\frac{1}{2}\\cos(2\\pi 3000 n t_s + \\pi/4) + \\cos(2\\pi 9000 n t_s)\n",
    "$$\n",
    "\n",
    "We can print each unique sample in $x(n)$ as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.around(a + b, 4)\n",
    "for i in range(0, N):\n",
    "    print(''.join(['x(', str(i), ') = ', str(x[i])]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Throughout this notebook, we will plot data using the `stem_plot` function below. This function will save a lot of notebook space and reduce code complexity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def stem_plot(x, y, title, xlabel, ylabel, xticks=None, yticks=None, figsize=(6, 3), subplot=(1, 1), bottom=None, style=None):\n",
    "    fig = plt.figure(figsize=figsize)\n",
    "    for idx, value in enumerate(y):\n",
    "        axes = fig.add_subplot(subplot[0], subplot[1], idx+1)\n",
    "        if bottom is not None:\n",
    "            axes.stem(x[idx], y[idx], use_line_collection = True, bottom=bottom[idx])\n",
    "        else:\n",
    "            axes.stem(x[idx], y[idx], use_line_collection = True)\n",
    "        if style is not None:\n",
    "            if style[idx] == 'dashed':\n",
    "                axes.plot(x[idx], y[idx], linestyle='dashed')\n",
    "        axes.grid(True, which='major')\n",
    "        if xticks is not None:\n",
    "            axes.set_xticks(xticks[idx])\n",
    "        if yticks is not None:\n",
    "            axes.set_yticks(yticks[idx])\n",
    "        axes.set_title(title[idx])\n",
    "        axes.set_xlabel(xlabel[idx])\n",
    "        axes.set_ylabel(ylabel[idx])\n",
    "        plt.box(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot $x(n)$ using the `stem_plot` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stem_plot(x      = [n],\n",
    "          y      = [x],\n",
    "          xticks = [n], \n",
    "          figsize= (6, 3),\n",
    "          title  = [\"Discrete Waveform $x(n)$\"], \n",
    "          style  = [\"dashed\"], \n",
    "          xlabel = [\"Samples (n)\"], \n",
    "          ylabel = [\"Amplitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. An Example of the DFT <a class=\"anchor\" id=\"an_example_of_the_dft\"></a>\n",
    "We can use the DFT to obtain the frequency content of a discrete waveform. The DFT equation is given below.\n",
    "\n",
    "$$\n",
    "X(k)=\\sum_{n=0}^{N-1}x(n)e^{-2\\pi kn/N}\n",
    "$$\n",
    "\n",
    "For a discrete input waveform $x(n)$ of length $N$, the DFT will output a discrete waveform $X(k)$ of the same length. Begin by initialising the output vector to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.zeros(N, dtype=np.complex64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compute the first output value $X(0)$, which corresponds to the DC offset or 0Hz component of $x(n)$.\n",
    "\n",
    "$$\n",
    "X(0) = \\sum_{n=0}^{15}x(n)e^{-2\\pi\\cdot 0 \\cdot n/16} = \\sum_{n=0}^{15}x(n)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[0] = np.sum(x)\n",
    "print(''.join(['X(0) = ', str(np.round(X[0], 2))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $x(n)$ does not have a DC offset, $X(0)$ is equal to zero. Now, calculate the second output value $X(1)$, which corresponds to a frequency component of 3kHz. We should expect this value to be non zero as our input waveform $x(n)$ contains a 3kHz cosine wave.\n",
    "\n",
    "$$\n",
    "X(1) = \\sum_{n=0}^{15}x(n)e^{-2\\pi \\cdot 1 \\cdot n/16}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[1] = np.sum(x*np.exp(-2j*np.pi*1*n/N))\n",
    "print(''.join(['X(1) = ', str(np.round(X[1], 2))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We continue by computing the third output value $X(2)$. This component corresponds to a frequency of 6kHz and should be equal to zero.\n",
    "\n",
    "$$\n",
    "X(2) = \\sum_{n=0}^{15}x(n)e^{-2\\pi \\cdot 2 \\cdot n/16}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[2] = np.sum(x*np.exp(-2j*np.pi*2*n/N))\n",
    "print(''.join(['X(2) = ', str(np.round(X[2], 2))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fourth output value $X(3)$ that corresponds to 9kHz should be more than zero as $x(n)$ contains a 9kHz cosine wave. This is computed as follows.\n",
    "$$\n",
    "X(3) = \\sum_{n=0}^{15}x(n)e^{-2\\pi \\cdot 3 \\cdot n/16}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[3] = np.sum(x*np.exp(-2j*np.pi*3*n/N))\n",
    "print(''.join(['X(3) = ', str(np.round(X[3], 2))]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than compute each value in $X(k)$, we can use a function to automate the process. We will name the function `dft` and use it throughout this notebook to compute the DFT as required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft(x):\n",
    "    N = np.size(x)\n",
    "    X = np.zeros(N, dtype=np.complex64)\n",
    "    n = np.arange(0, N)\n",
    "    for k in range(N):\n",
    "        X[k] = np.sum(x*np.exp(-2j*np.pi*k*n/N))\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, compute the DFT of $x(n)$ and print the results for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.round(dft(x), 2) + (0.0 + 0.0j) # Add 0.0+0.0j to remove -0.0-0.0j\n",
    "for i in range(N):\n",
    "    print(''.join(['X(',str(i),') = ',str(X[i])]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown, there appears to be symmetry in the results printed above. For instance, $|X(1)| = |X(15)|$ and $|X(3)| = |X(13)|$. We can plot the results in the following sections to aid our understanding of the DFT output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Magnitude and Phase <a class=\"anchor\" id=\"magnitude_and_phase\"></a>\n",
    "We can obtain the real and imaginary components of the DFT output by using the equations below.\n",
    "\n",
    "$$\n",
    "\\text{Re}(c) = \\frac{c + \\overline{c}}{2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{Im}(c) = \\frac{c - \\overline{c}}{2j}\n",
    "$$\n",
    "\n",
    "The magnitude of the DFT output can then be obtained using,\n",
    "\n",
    "$$\n",
    "|X(k)| = \\sqrt{\\text{Re}(X(k))^{2}+\\text{Im}(X(k))^2}.\n",
    "$$\n",
    "\n",
    "We can simply compute the magnitude of a complex vector using `np.abs`. Then, round the vector to two decimal places."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xm = np.round(np.abs(X), 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will plot the magnitude in a moment. We can also obtain the phase of the DFT output as,\n",
    "\n",
    "$$\n",
    "X_{\\text{Ø}}(k) = \\tan^{-1}\\left(\\frac{\\text{Im}(X(k))}{\\text{Re}(X(k))}\\right).\n",
    "$$\n",
    "\n",
    "This calculation can be performed using `np.angle` (also rounding the output to two decimal places)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xp = np.round(np.angle(X, deg=True), 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DFT output obtained in the previous section can now be represented in terms of its magnitude and phase. These values can be plotted for inspection. Again, you will notice symmetry in the output result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in range(0, N):\n",
    "    print(''.join(['X(',str(i),') = ',str(Xm[i]),' ∠ ',str(Xp[i]),'°']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Plotting the DFT <a class=\"anchor\" id=\"plotting_the_dft\"></a>\n",
    "Plotting the DFT output is a very useful when analysing the frequency content of a discrete waveform. The real and imaginary components of $X(k)$ can be plotted as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stem_plot(x      = [n*(fs/N)/1e3, n*(fs/N)/1e3],\n",
    "          y      = [np.real(X), np.imag(X)],\n",
    "          yticks = [np.arange(0, 9), np.arange(-3, 4)],\n",
    "          figsize= (14, 3),\n",
    "          xticks = [n*(fs/N)/1e3, n*(fs/N)/1e3], \n",
    "          subplot= (1, 2),\n",
    "          title  = [\"Real Part of $X(k)$\", \"Imaginary Part of $X(k)$\"],\n",
    "          xlabel = [\"Frequency (kHz)\", \"Frequency (kHz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Amplitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Symmetry is apparent around the y-axis in the plot showing the real part of $X(k)$ (excluding the 0kHz component). Similarly, we can also see that the imaginary part of $X(k)$ is almost symmetric about the 24kHz (again ignoring the 0KHz component).\n",
    "\n",
    "We can also plot the magnitude and phase of $X(k)$ that we computed in the previous section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stem_plot(x      = [n*(fs/N)/1e3, n*(fs/N)/1e3],\n",
    "          y      = [Xm, Xp],\n",
    "          yticks = [np.arange(0, 9), np.arange(-45, 60, 15)],\n",
    "          figsize= (14, 3),\n",
    "          xticks = [n*(fs/N)/1e3, n*(fs/N)/1e3], \n",
    "          subplot= (1, 2),\n",
    "          title  = [\"Magnitude of $X(k)$\", \"Phase of $X(k)$\"],\n",
    "          xlabel = [\"Frequency (kHz)\", \"Frequency (kHz)\"], \n",
    "          ylabel = [\"Magnitude\", \"Phase (degrees)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can normalise the output magnitude of $X(k)$ by dividing it by $N$. Note, that this technique only works for rectangular windows (discussed later).\n",
    "\n",
    "$$\n",
    "X'(k) = \\frac{1}{N}\\sum_{n=0}^{N-1}x(n)e^{-j2\\pi kn/N}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xnorm = (Xm/N)\n",
    "\n",
    "stem_plot(x      = [n*(fs/N)/1e3],\n",
    "          y      = [Xnorm],\n",
    "          figsize= (6, 3),\n",
    "          yticks = [np.arange(0, 1.125, 0.125)],\n",
    "          xticks = [n*(fs/N)/1e3],\n",
    "          title  = [\"Normalised Magnitude of $X(k)$\"],\n",
    "          xlabel = [\"Frequency (kHz)\"], \n",
    "          ylabel = [\"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the normalised magnitude frequency spectrum above contains lobes at 3kHz and 9kHz, which corresponds with the original cosine waveforms in $x(n)$. However, the value of the 3kHz lobe is 0.25 and the 9kHz lobe is 0.5. This problem can be resolved by considering positive and negative frequencies, and the periodicity of the DFT in the following section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Positive and Negative Frequencies <a class=\"anchor\" id=\"symmetry_and_periodicity\"></a>\n",
    "All real waveforms have a positive and negative complex exponential of the same frequency. For example, we can use Euler's formula to describe $x(n)$ using complex exponentials.\n",
    "\n",
    "$$\n",
    "x(n) = \\frac{1}{2}\\left(\\frac{e^{j2\\pi 3000 n t_s + \\pi/4} + e^{-j2\\pi 3000 n t_s + \\pi/4}}{2}\\right) + \\left(\\frac{e^{j2\\pi 9000 n t_s} + e^{-j2\\pi 9000 n t_s}}{2}\\right)\n",
    "$$\n",
    "\n",
    "Rearranging the equation above so that positive and negative complex exponentials are paired together gives,\n",
    "\n",
    "$$\n",
    "x(n) = \\frac{1}{4}e^{j2\\pi 3000 n t_s + \\pi/4} + \\frac{1}{2}e^{j2\\pi 9000 n t_s} + \\frac{1}{4}e^{-j2\\pi 3000 n t_s + \\pi/4} + \\frac{1}{2}e^{-j2\\pi 9000 n t_s}.\n",
    "$$\n",
    "\n",
    "Plotting the frequencies of the complex exponentials against their amplitude produces the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stem_plot(x      = [(n-N/2+1)*(fs/N)/1e3],\n",
    "          y      = [np.concatenate((Xnorm[N//2+1:N], Xnorm[0:N//2+1]))],\n",
    "          figsize= (6, 3),\n",
    "          yticks = [np.arange(0, 1.25, 0.25)],\n",
    "          xticks = [(n-N/2+1)*(fs/N)/1e3],\n",
    "          title  = [\"Complex Frequency Spectrum of $|X\\'(k)|$\"],\n",
    "          xlabel = [\"Frequency (kHz)\"], \n",
    "          ylabel = [\"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the plot above is similar to the normalised magnitude of $X(k)$ in the previous section. However, this plot shows positive and negative frequency components. We can sum the positive and negative components together to obtain a real waveform. For instance, let us demonstrate this process using the 9kHz cosine wave from before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "positive_exp = (1/2) * np.exp(2j*np.pi*9000*n*ts)\n",
    "negative_exp = (1/2) * np.exp(-2j*np.pi*9000*n*ts)\n",
    "real_waveform = positive_exp + negative_exp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now present the complex exponentials and their output sum using 3-dimensional complex plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "fig = plt.figure(figsize=(14,5))\n",
    "title = ['Positive-Frequency\\nComplex Exponential',\n",
    "         'Negative-Frequency\\nComplex Exponential',\n",
    "         'Sum of Complex\\nExponentials (Cosine Wave)']\n",
    "\n",
    "for idx, value in enumerate([positive_exp, negative_exp, real_waveform]):\n",
    "    axes = fig.add_subplot(1,3,idx+1, projection='3d')\n",
    "    axes.plot(np.real(value), n, np.imag(value))\n",
    "    axes.set_proj_type('ortho')\n",
    "    axes.set_xticks(np.arange(-1, 1.5, 0.5))\n",
    "    axes.set_zticks(np.arange(-1, 1.5, 0.5))\n",
    "    axes.set_yticks(np.arange(0, N+1, N/4))\n",
    "    axes.set_xlabel('Re')\n",
    "    axes.set_ylabel('Samples (n)')\n",
    "    axes.set_zlabel('Im')\n",
    "    axes.set_title(title[idx])\n",
    "    axes.view_init(-150, 45)\n",
    "    axes.dist = 12\n",
    "    \n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots above show the positive-frequency complex exponential rotating in one direction, and the negative-frequency complex exponential rotating in the other direction. The sum of these complex exponential produces the original 9kHz real valued waveform.\n",
    "\n",
    "Negative frequencies are possible if we consider the periodicity of the DFT using the equation below.\n",
    "\n",
    "$$\n",
    "X(k) = X(k+N) = \\sum_{n=0}^{N-1}x(n)e^{-j2\\pi kn/N}\n",
    "$$\n",
    "\n",
    "We can compute the DFT when $k$ operates across the range $[-16, 15]$ and compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(N):\n",
    "    temp = np.round(np.sum(x*np.exp(-2j*np.pi*i*n/N)), 2) + (0.0 + 0.0j)\n",
    "    temp_N = np.round(np.sum(x*np.exp(-2j*np.pi*(i-N)*n/N)), 2) + (0.0 + 0.0j)\n",
    "    print('{:s}{:d}{:2s}{:2s}{:.2f}{:>4s}{:d}{:2s}{:2s}{:.2f}'.format(\n",
    "        'X(', i-N, ')', '=', temp_N, 'X(', i, ')', '=', temp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the DFT output is periodic every $N$ samples. We normally multiply the output DFT by 2 for real waveforms only and can simply ignore one half of the DFT as it contains redundant data.\n",
    "\n",
    "It is important to note that complex waveforms do not have a combination of positive and negative complex exponentials. They are characterised by one complex exponential only. Consider the complex wave given by\n",
    "\n",
    "$$\n",
    "x(n) = \\frac{1}{4}e^{-j2\\pi 12000 n t_s}.\n",
    "$$\n",
    "\n",
    "We can compute the DFT and plot the results for analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = -12e3 # select from -21e3 to 24e3 in multiples of 3e3\n",
    "cwave = np.round((1/4)*np.exp(2j*np.pi*fc*n*ts), 2)\n",
    "cwave_fnorm = np.round(np.abs(dft(cwave)), 2)/N\n",
    "fmag_shift = np.concatenate((cwave_fnorm[N//2+1:N], cwave_fnorm[0:N//2+1]))\n",
    "\n",
    "stem_plot(x      = [(n-N/2+1)*(fs/N)/1e3],\n",
    "          y      = [fmag_shift],\n",
    "          figsize= (6, 3),\n",
    "          yticks = [np.arange(0, 0.625, 0.125)],\n",
    "          xticks = [(n-N/2+1)*(fs/N)/1e3],\n",
    "          title  = [\"Complex Frequency Spectrum\"],\n",
    "          xlabel = [\"Frequency (kHz)\"], \n",
    "          ylabel = [\"Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this complex waveform only contains one frequency component in its DFT output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Magnitude, Power, and the Logarithmic Scale <a class=\"anchor\" id=\"magnitude_power_log\"></a>\n",
    "The DFT output can be plotted using several different units and scales. Consider the discrete time domain waveform given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1000\n",
    "ts = 1/fs\n",
    "fd = 200\n",
    "N = 32\n",
    "qlevels = 128 # -128 to 127\n",
    "n = np.arange(0, N)\n",
    "x = qlevels * np.sin(2*np.pi*fd*ts*n)\n",
    "\n",
    "stem_plot(x      = [n],\n",
    "          y      = [x],\n",
    "          xticks = [n], \n",
    "          yticks = [np.arange(-qlevels, qlevels+1, qlevels//2)],\n",
    "          figsize= (14, 3),\n",
    "          title  = [\"Quantised Discrete Waveform $x(n)$\"], \n",
    "          style  = [\"dashed\"], \n",
    "          xlabel = [\"Samples (n)\"], \n",
    "          ylabel = [\"Amplitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power spectrum of this discrete waveform can be obtained using,\n",
    "\n",
    "$$\n",
    "X'_{\\text{ps}}(k)=|X'(k)|^2.\n",
    "$$\n",
    "\n",
    "Similarly, we can also represent the power spectrum using a log-scale equivalent, as given below.\n",
    "\n",
    "$$\n",
    "X'_{\\text{dB}}(k)=20\\text{log}_{10}(|X'(k)|) \\, \\text{dB}\n",
    "$$\n",
    "\n",
    "We can plot the power spectrum using the following code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.round(dft(x), 2) + (0.0 + 0.0j) # Add 0.0+0.0j to remove -0.0-0.0j\n",
    "Xm = np.round(np.abs(X), 2)\n",
    "Xnorm = (Xm/N)*2 # multiply by 2 as real\n",
    "Xpow = Xnorm**2  # power calculation\n",
    "Xlogpow = 20*np.where(Xnorm>0, np.log10(Xnorm), 0) # log scale\n",
    "\n",
    "stem_plot(x      = [n[0:N//2+1:1]*(fs/N), n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [Xpow[0:N//2+1:1], Xlogpow[0:N//2+1:1]],\n",
    "          subplot= (1, 2),\n",
    "          figsize= (14, 3),\n",
    "          bottom = [0, 15],\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(0, 10000+1, 2000), np.arange(15, 40+1, 5)],\n",
    "          xticks = [n[0:N//2+1:2]*fs/N, n[0:N//2+1:2]*fs/N],\n",
    "          title  = [\"Power Spectra of $X(k)$\", \"Log-Scale Power Spectra of $X(k)$\"],\n",
    "          xlabel = [\"Frequency (Hz)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Power\", \"Power (dB)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to use dBFS (decibels relative to full-scale) to represent the power spectrum, which is given as,\n",
    "\n",
    "$$\n",
    "X'_{\\text{dBFS}}(k)=20\\text{log}_{10}\\left(\\frac{|X'(k)|}{\\text{full-scale}}\\right) \\, \\text{dBFS},\n",
    "$$\n",
    "\n",
    "where full-scale is the maximum quantisation level of a discrete waveform. We can plot the dBFS representation of the discrete time waveform below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xdbfs = 20*np.where(Xpow>0, np.log10(Xnorm/qlevels), 0)\n",
    "\n",
    "stem_plot(x      = [n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [Xdbfs[0:N//2+1:1]],\n",
    "          figsize= (6, 3),\n",
    "          bottom = [-25],\n",
    "          yticks = [np.arange(0, -25-1, -5)],\n",
    "          style  = [\"dashed\"],\n",
    "          xticks = [n[0:N//2+1:2]*fs/N],\n",
    "          title  = [\"dBFS Log-Scale Power Spectra of $X(k)$\"],\n",
    "          xlabel = [\"Frequency (Hz)\"],\n",
    "          ylabel = [\"Power (dBFS)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Frequency Bins and Resolution <a class=\"anchor\" id=\"frequency_bins\"></a>\n",
    "The frequency axis of the above spectrum plots are divided into several discrete frequencies, which are known as bins. The spacing between the bins is known as the frequency resolution. We can use an example to demonstrate how the DFT length $N$, and the sample frequency $f_s$ affect the frequency resolution of the DFT output.\n",
    "\n",
    "Consider the discrete sine wave given below, which has a sample frequency of 1000Hz and contains 8 samples. The target frequency of the sine wave is 125Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1000\n",
    "ts = 1/fs\n",
    "fd = 125\n",
    "N = 8\n",
    "n = np.arange(0, N)\n",
    "x = np.sin(2*np.pi*fd*ts*n)\n",
    "X = dft(x)\n",
    "Xm = np.round(np.abs(X), 2)\n",
    "Xnorm = (Xm/N)*2 # multiply by 2 as real\n",
    "\n",
    "stem_plot(x      = [n, n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [x, Xnorm[0:N//2+1:1]],\n",
    "          subplot= (1, 2),\n",
    "          figsize= (14, 3),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 1.2, 0.2)],\n",
    "          xticks = [n, n[0:N//2+1:1]*fs/N],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Normalised Magnitude of $X(k)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown above, the normalised magnitude spectrum consists of 5 evenly spaced frequencies, which are 0Hz, 125Hz, 250Hz, 375Hz, and 500Hz. We can see that the frequency resolution is 125Hz. The frequency resolution, denoted as $\\Delta f$ can also be calculated using,\n",
    "\n",
    "$$\n",
    "\\Delta f = \\frac{f_s}{N} = \\frac{1000 \\text{ samples/second}}{8 \\text{ samples}} = 125 \\text{ Hz}.\n",
    "$$\n",
    "\n",
    "The frequency resolution is affected by the number of samples in the input waveform. For instance, increasing the number of samples such that $N=16$ improves the frequency resolution, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 16\n",
    "n = np.arange(0, N)\n",
    "x = np.sin(2*np.pi*fd*ts*n)\n",
    "X = dft(x)\n",
    "Xm = np.round(np.abs(X), 2)\n",
    "Xnorm = (Xm/N)*2 # multiply by 2 as real\n",
    "\n",
    "stem_plot(x      = [n, n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [x, Xnorm[0:N//2+1:1]],\n",
    "          subplot= (1, 2),\n",
    "          figsize= (14, 3),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 1.2, 0.2)],\n",
    "          xticks = [n, n[0:N//2+1:1]*fs/N],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Normalised Magnitude of $X(k)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the frequency resolution improves if the sample frequency is reduced. For example, the code cell below uses a sample frequency of $f_s=500\\text{Hz}$ and $N=8$. Notice that the frequency resolution is now $\\Delta f = 62.5\\text{Hz}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 500\n",
    "ts = 1/fs\n",
    "N = 8\n",
    "n = np.arange(0, N)\n",
    "x = np.sin(2*np.pi*fd*ts*n)\n",
    "X = dft(x)\n",
    "Xm = np.round(np.abs(X), 2)\n",
    "Xnorm = (Xm/N)*2 # multiply by 2 as real\n",
    "\n",
    "stem_plot(x      = [n, n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [x, Xnorm[0:N//2+1:1]],\n",
    "          subplot= (1, 2),\n",
    "          figsize= (14, 3),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 1.2, 0.2)],\n",
    "          xticks = [n, n[0:N//2+1:1]*fs/N],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Normalised Magnitude of $X(k)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can improve the frequency resolution by decreasing the sample rate and/or increasing the number of samples used by the DFT. However, we must consider the implications of changing these parameters. Increasing the number of samples $N$ means that additional computation is necessary to compute the DFT. Also, decreasing the sample frequency has the effect of reducing the range of frequencies that can be observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Spectral Leakage <a class=\"anchor\" id=\"spectral_leakage\"></a>\n",
    "The DFT operates very well on discrete waveforms that contain frequency components that are an integer multiple of the frequency resolution (or fundamental frequency). When the input waveform does not fall directly into a DFT bin, spectral leakage occurs. For instance, consider an 80Hz discrete wave that has been sampled at 1000Hz and contains 16 samples. The corresponding DFT does not have a frequency bin for the 80Hz component, therefore, the energy spreads into neighbouring frequencies. This phenomenon is known to cause many people to interpret DFT results incorrectly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1000\n",
    "ts = 1/fs\n",
    "fd = 80\n",
    "N = 16\n",
    "n = np.arange(0, N)\n",
    "x = np.sin(2*np.pi*fd*ts*n)\n",
    "X = dft(x)\n",
    "Xm = np.round(np.abs(X), 2)\n",
    "Xnorm = (Xm/N)*2 # multiply by 2 as real\n",
    "\n",
    "stem_plot(x      = [n, n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [x, Xnorm[0:N//2+1:1]],\n",
    "          subplot= (1, 2),\n",
    "          figsize= (14, 3),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 1.2, 0.2)],\n",
    "          xticks = [n, n[0:N//2+1:1]*fs/N],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Normalised Magnitude of $X(k)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implicitly assume that the input waveform to the DFT contains pseudo-periods every $N$ samples. In the discrete waveform above, there will be an abrupt transition between one period and another. These abrupt changes in amplitude are known as discontinuities and are only found at the endpoints of a waveform's sampling interval.\n",
    "\n",
    "It is possible to plot the discrete magnitude response we obtained above, alongside an approximation of the continuous frequency response of the 80Hz sine wave. This can be seen by running the code cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nc = 1024\n",
    "Xc = dft(np.pad(x, (0, Nc-N), 'constant'))\n",
    "Xcm = np.round(np.abs(Xc), 2)\n",
    "Xcnorm = (Xcm/N)*2 # multiply by 2 as real\n",
    "\n",
    "fig = plt.figure(figsize=(6, 3))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.stem(n[0:N//2+1:1]*(fs/N), Xnorm[0:N//2+1:1], use_line_collection=True)\n",
    "axes.plot(np.arange(0,Nc//2+1,1)*(fs/Nc), Xcnorm[0:Nc//2+1:1], linestyle='dashed')\n",
    "axes.grid(True, which='major')\n",
    "axes.set_title('Normalised Magnitude of $X(k)$')\n",
    "axes.set_yticks(np.arange(0, 1.2, 0.2))\n",
    "axes.set_xticks(n[0:N//2+1:1]*fs/N)\n",
    "axes.set_xlabel('Frequency (kHz)')\n",
    "axes.set_ylabel('Normalised Magnitude')\n",
    "plt.box(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main lobe of the continuous frequency response is centred across 80Hz. As shown, the discrete response is simply a sampled version of the continuous frequency magnitude response. Each bin of the discrete response has been given a value that is an equivalent point on the continuous magnitude response."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Windowing <a class=\"anchor\" id=\"windowing\"></a>\n",
    "Spectral leakage can be suppressed by windowing a discrete waveform before using the DFT. Windowing is the process of extracting a subset of data from a larger set (extracting an interval of data). Many windows in Digital Signal Processing (DSP) have the effect of \"tapering\" the endpoints of a discrete waveform to reduce the possibility of abrupt transitions, which cause spectral leakage to occur. We can describe windowing in the DFT equation as below, where $w(n)$ denotes the discrete window.\n",
    "\n",
    "$$\n",
    "X_w(k)=\\sum_{n=0}^{N-1}w(n) \\cdot x(n)e^{-j2\\pi kn/N}\n",
    "$$\n",
    "\n",
    "In this section, we will briefly exploring windowing. We will begin by creating a test waveform for analysis. The sine wave below has a sampling frequency of 1000Hz, a desired frequency of 80Hz and a sample length of 16."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1000\n",
    "ts = 1/fs\n",
    "fd = 80\n",
    "N = 16\n",
    "n = np.arange(0, N)\n",
    "x = np.sin(2*np.pi*fd*ts*n)\n",
    "\n",
    "stem_plot(x      = [n],\n",
    "          y      = [x],\n",
    "          figsize= (6, 3),\n",
    "          xticks = [n], \n",
    "          yticks = [np.arange(-1, 1.5, 0.5)],\n",
    "          title  = [\"Discrete Waveform $x(n)$\"], \n",
    "          style  = [\"dashed\"], \n",
    "          xlabel = [\"Samples (n)\"], \n",
    "          ylabel = [\"Amplitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest form of window is known as a rectangular window, which does not modify (or taper) the extracted data. Rectangular windows are used to extract an interval of data from a discrete waveform. Since the waveform above has not been modified, it has been extracted using a rectangular window of length 16.\n",
    "\n",
    "A useful window for tapering the endpoints of a discrete waveform is known as a Hanning window. You can see its effects when applied to the data above by running the following code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.hanning(N)\n",
    "xw = x*w\n",
    "\n",
    "stem_plot(x      = [n],\n",
    "          y      = [xw],\n",
    "          figsize= (6, 3),\n",
    "          xticks = [n], \n",
    "          yticks = [np.arange(-1, 1.5, 0.5)],\n",
    "          title  = [\"Windowed Discrete Waveform $x(n)$\"], \n",
    "          style  = [\"dashed\"], \n",
    "          xlabel = [\"Samples (n)\"], \n",
    "          ylabel = [\"Amplitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the effects of using a rectangular window and Hanning window by observing their corresponding magnitude spectra plots, which are generated using the DFT. Run the code cell below to generate these plots. You will notice that the rectangular window will have a larger amount of spectral leakage than the Hanning window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xm = np.round(np.abs(dft(x)), 2)*2\n",
    "Xwm = np.round(np.abs(dft(xw)), 2)*2\n",
    "\n",
    "stem_plot(x      = [n[0:N//2+1:1]*(fs/N), n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [Xm[0:N//2+1:1], Xwm[0:N//2+1:1]],\n",
    "          yticks = [np.arange(0, 20, 4), np.arange(0, 20, 4)],\n",
    "          figsize= (14, 3),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          xticks = [n[0:N//2+1:1]*fs/N, n[0:N//2+1:1]*fs/N], \n",
    "          subplot= (1, 2),\n",
    "          title  = [\"Magnitude Spectra (Rectangular Window)\",\n",
    "                    \"Magnitude Spectra (Hanning Window)\"],\n",
    "          xlabel = [\"Frequency (Hz)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Magnitude\", \"Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the Hanning window exhibits less spectral leakage than the rectangular window, it has a wider main lobe. We can normalise the effects of windows by dividing the output DFT by the sum of the window, which is given as\n",
    "\n",
    "$$\n",
    "X_w'(k)=\\frac{1}{\\sum_{n}{w(n)}}X_w(k).\n",
    "$$\n",
    "\n",
    "Normalising each DFT output by their corresponding window produces the following normalised magnitude spectra plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xnorm = Xm/np.sum(np.ones(N))\n",
    "Xwnorm = Xwm/np.sum(np.hanning(N))\n",
    "\n",
    "stem_plot(x      = [n[0:N//2+1:1]*(fs/N), n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [Xnorm[0:N//2+1:1], Xwnorm[0:N//2+1:1]],\n",
    "          yticks = [np.arange(0, 1.25, 0.25), np.arange(0, 1.25, 0.25)],\n",
    "          figsize= (14, 3),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          xticks = [n[0:N//2+1:1]*fs/N, n[0:N//2+1:1]*fs/N], \n",
    "          subplot= (1, 2),\n",
    "          title  = [\"Normalised Magnitude Spectra (Rectangular Window)\",\n",
    "                    \"Normalised Magnitude Spectra (Hanning Window)\"],\n",
    "          xlabel = [\"Frequency (Hz)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Normalised Magnitude\", \"Normalised Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is clear that tapering the endpoints of a discrete waveform before it undergoes the DFT is useful to suppress spectral leakage. Using a non-rectangular window to taper a discrete waveform's endpoints always has the effect of increasing the width of the primary lobe.\n",
    "\n",
    "We can plot many different windows from the Numpy Python library using the code cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 16 # even number only\n",
    "windows = ['ones', 'hamming', 'hanning', 'blackman', 'bartlett']\n",
    "names = ['Rectangular', 'Hamming', 'Hanning', 'Blackman', 'Bartlett']\n",
    "\n",
    "fig = plt.figure(figsize=(14, 16))\n",
    "\n",
    "for idx, window in enumerate(windows):\n",
    "    w = getattr(np, window)\n",
    "    axes = fig.add_subplot(3, 2, idx+1)\n",
    "    if window == 'bartlett':\n",
    "        n = np.arange(0, N+1, 1)\n",
    "        axes.stem(n, w(N+1), use_line_collection=True)\n",
    "        axes.plot(n, w(N+1), linestyle='dashed')\n",
    "    else:\n",
    "        n = np.arange(0, N, 1)\n",
    "        axes.stem(n, w(N), use_line_collection=True)\n",
    "        axes.plot(n, w(N), linestyle='dashed')\n",
    "    axes.grid(True, which='major')\n",
    "    axes.set_xticks(n)\n",
    "    axes.set_title(''.join([names[idx], ' Window']))\n",
    "    axes.set_xlabel('Samples (n)')\n",
    "    axes.set_ylabel('Amplitude')\n",
    "    plt.box(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each window has their own magnitude response in the frequency domain. These can be observed by running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning) \n",
    "\n",
    "N = 64\n",
    "Nc = 16384\n",
    "fig = plt.figure(figsize=(14, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "colours = [\"b\", \"r\", \"y\", \"m\", \"g\"]\n",
    "\n",
    "for idx, window in enumerate(windows):\n",
    "    w = getattr(np, window)\n",
    "    Xw = np.fft.fft(np.pad(w(N), (0, Nc-N), 'constant'))\n",
    "    Xwm = np.abs(Xw)\n",
    "    Xwnorm = Xwm/np.sum(w(N))\n",
    "    Xwlog = 20*np.log10(Xwnorm)\n",
    "    axes.plot(np.arange(0, Nc//8+1, 1)/8192, Xwlog[0:Nc//8+1:1], colours[idx])\n",
    "    axes.set_yticks(np.arange(0, -90, -10))\n",
    "    axes.set_xticks(np.arange(0, N/(Nc/N)+0.05, 0.05))\n",
    "    axes.set_ylim([-80, 0])\n",
    "    axes.grid(True, which='major')\n",
    "    axes.legend(names)\n",
    "    plt.box(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the width of the rectangular window's primary lobe is the smallest of all windows above. However, the subsequent side lobes have very little suppression. The other tapered windows share a common trait of a wide primary lobe and suppressed side lobes (of different levels)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. Scalloping Loss <a class=\"anchor\" id=\"scalloping_loss\"></a>\n",
    "We can measure the worst-case reduction in energy when a discrete waveform contains a frequency component that resides precisely between two DFT bins. This measurement is known as Scalloping Loss (SL) [1] and is defined as\n",
    "\n",
    "$$\n",
    "\\text{SL} = \\frac{|\\sum_{n}w(n)e^{-j\\pi n/N}|}{\\sum_{n}w(n)}.\n",
    "$$\n",
    "\n",
    "We can define a simple Python function that can measure the SL of any discrete window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scallop_loss(w):\n",
    "    N = len(w)\n",
    "    n = np.arange(0, N, 1)\n",
    "    numerator = np.abs(np.sum(w*np.exp(-1j*np.pi*n/N)))\n",
    "    decimator = np.sum(w)\n",
    "    return numerator/decimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can use the function to measure the SL of a rectangular window of length 16."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(scallop_loss(np.ones(16)), 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the SL can be measured for a Hanning window of the same length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(scallop_loss(np.hanning(16)), 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the SL for the rectangular window is worse than the SL for the Hanning window. The Hanning window has a wider main lobe than a rectangular window, so these results are expected.\n",
    "\n",
    "You can compute the SL for any window using the function above. As an exercise, you can compute the SL for all windows, and determine the window that exhibits the best SL."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Zero Padding <a class=\"anchor\" id=\"zero_padding\"></a>\n",
    "The final topic in this notebook is the zero padding technique, which is very useful when improving the frequency resolution of the DFT output. Consider the discrete waveform and its corresponding DFT output given in the cell below. The discrete wave was sampled using 2000Hz, the desired frequency is 250Hz, and it has 16 samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 2000\n",
    "ts = 1/fs\n",
    "fd = 250\n",
    "N = 16\n",
    "n = np.arange(0, N, 1)\n",
    "x = np.sin(2*np.pi*fd*ts*n)\n",
    "Xm = np.round(np.abs(dft(x)), 2)*2\n",
    "\n",
    "stem_plot(x      = [n, n[0:N//2+1:1]*(fs/N)],\n",
    "          y      = [x, Xm[0:N//2+1]],\n",
    "          subplot= (2, 1),\n",
    "          figsize= (14, 10),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 20, 4)],\n",
    "          xticks = [n, n[0:N//2+1:1]*fs/N],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Magnitude Spectra of $x(n)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frequency resolution of the DFT above is $\\Delta f = 125\\text{Hz}$. We can zero pad the end of the discrete wave with 16 zeros and then perform a 32-point DFT to improve the frequency resolution. This is implemented in the code cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nc = 32\n",
    "n = np.arange(0, Nc, 1)\n",
    "xp = np.pad(x, (0, Nc-N), 'constant')\n",
    "Xm = np.round(np.abs(dft(xp)), 2)*2\n",
    "\n",
    "stem_plot(x      = [n, n[0:Nc//2+1:1]*(fs/Nc)],\n",
    "          y      = [xp, Xm[0:Nc//2+1]],\n",
    "          subplot= (2, 1),\n",
    "          figsize= (14, 10),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 20, 4)],\n",
    "          xticks = [n, n[0:Nc//2+1:1]*fs/Nc],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Magnitude Spectra of $x(n)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the frequency resolution has increased to $\\Delta f = 62.5\\text{Hz}$. We can append another 16 zeros to the end of the signal to improve the frequency resolution further, as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nc = 48\n",
    "n = np.arange(0, Nc, 1)\n",
    "xp = np.pad(x, (0, Nc-N), 'constant')\n",
    "Xm = np.round(np.abs(dft(xp)), 2)*2\n",
    "\n",
    "stem_plot(x      = [n, n[0:Nc//2+1:1]*(fs/Nc)],\n",
    "          y      = [xp, Xm[0:Nc//2+1]],\n",
    "          subplot= (2, 1),\n",
    "          figsize= (14, 10),\n",
    "          style  = [\"dashed\", \"dashed\"],\n",
    "          yticks = [np.arange(-1, 1+0.5, 0.5), np.arange(0, 20, 4)],\n",
    "          xticks = [n, n[0:Nc//2+1:1]*fs/Nc],\n",
    "          title  = [\"Discrete Waveform $x(n)$\", \"Magnitude Spectra of $x(n)$\"],\n",
    "          xlabel = [\"Samples (n)\", \"Frequency (Hz)\"], \n",
    "          ylabel = [\"Amplitude\", \"Magnitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the output magnitude spectra now has a frequency resolution of $\\Delta f = 41.67\\text{Hz}$.\n",
    "\n",
    "Keep in mind that zero padding has the effect of improving the frequency resolution of the DFT, but increases its computational requirements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Conclusion <a class=\"anchor\" id=\"conclusion\"></a>\n",
    "This notebook has presented the DFT and demonstrated several examples showcasing its underlying operation. We investigated DFT topics including periodicity, positive and negative frequencies, frequency representation, spectral leakage, windowing, scalloping loss, and zero padding.\n",
    "\n",
    "In the next notebook, we introduce the Fast Fourier Transform (FFT) and spectrograms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](03_fouriers_theorem.ipynb) || [Next Notebook 🚀](05_fast_fourier_transform.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
